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The solution of the Ornstein-Zernike equation with various closure ap- 
proximations is studied. This problem is rewritten as an integral equation 
that can be solved iteratively on a grid. The convergence of the fixed point 
iterations is relatively slow. We consider transformations of the sequence 
of solution vectors using non-linear sequence transformations, so-called 
vector extrapolation processes. An example is the vector J transforma- 
tion. The transformed vector sequences turn out to converge considerably 
faster than the original sequences. 



1 Classical Many-Particle Systems 

In this paper we investigate acceleration methods for solving the fundamen- 
tal equation for the pair distribution function of classical many-particle sys- 
tems, the so-called Ornstein-Zernike equation. The thermodynamic proper- 
ties of such systems are determined by the interaction between the particles 
from which the system is built up. If one knows the two-particle distribution 
function </, one can calculate all thermodynamic properties of the considered 
system, g is defined in the canonical ensemble by [1, Chapter 4] 
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where V is the volume of the system, n is the number of particles, and as 
usual /3 = [kBT)~^^ where is Boltzmann's constant and T is the absolute 
temperature. Note that the denominator is the classical configuration integral. 
We restrict our attention to pair potentials u, i.e., 

n 

[/(ri,...,r„) = ^M(r„rj) . (2) 

i<j 



For simplicity we consider in a first step only systems with radially symmetric 
interactions between identical particles. For the theoretical development (see 
e.g. [1, Chapters 6, 7]) of the equations it is useful to define the Mayer /- 
function: 

/(r):=e-'^"M-l, (3) 



where u[r) is the potential energy between particle 1 and 2 at distance r. The 
latter is said to be regular (short ranged) (see [2, p. 72]) if it is bounded below 
and satisfies 



dT{r) < oo V /3 > 0. (4) 



Of this type are for example the Lennard- Jones (LJ) potential 



where cr is a distance parameter and e is the depth of the potential, and the 
hard sphere potential 

oo V r < cr 

u{r) = <^ . (6) 

V r > (7 



On the other hand, there are pair potentials u[r) which do not obey rela- 
tion (4). Nevertheless, they lead to thermodynamical behavior of systems of 
particles interacting with such u[r). A famous example is the classical one- 
component plasma (OCP) with the pair potential 

<-) = ^^k^T-^ r, = ^^^^-^-^, -=(^1 (7) 



for particles of charge Z in a neutralizing background. Here, eo is the absolute 
value of the elementary charge, eo is the dielectric constant, the relative 
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permittivity, and p is the average number density that can be used to define 
a length scale a. The plasma parameter Vp is a dimensionless quantity. For 
further convenience we divide the potential in a long-range part u^^\r) and a 
short-range part u^^\r) in the following manner according to [3]: 

u^\r) ■=V.pkBT-eii{ar) , u^'\r) := u{r) - u^\r) , (8) 
r 

where a is a parameter to be chosen (usually a = 1.08/a, see [3]). For the 
definition of the error function erf(a;) see [4, Chapter 7]. The Fourier transform 
of u^^\r) can - similarly to the Fourier transform of the Coulomb potential - 
be calculated in the distributional sense. It is short ranged and given by 

= 47rr,A;Br^ exp(--^). (9) 

If the potential u is radially symmetric and therefore only a function of r : = 
|ri — we can establish the pair distribution function function of r. 

To determine this quantity we are using the Ornstein-Zernike (OZ) equation 

[1]: 

h = c + p-c*h (10) 

where * denotes a convolution defined by 

[/*5](r) = / f{v-v')g{v')dT{v'). (11) 

The density p is the average number density. The function h{r) := g{r) — 1 is 
called the total correlation function and c(r) the direct correlation function. 
We note that the convolution of two radially symmetric functions is again a 
radially symmetric function. For the two unknown functions h and c we need 
a second equation, which is called the closure of the OZ equation and is given 
in general by [1]: 

g{r) = exp(— /3M(r) -|- h{r) — c[r) -\- E{r)) . (12) 

E is an infinite sum of multicenter integrals, the so called bridge diagrams, 
which are known in principle as complicated multidimensional integrals. These 
are very hard to evaluate. Thus, usually various simple approximations are 
used for them. E{r) = is the HyperNetted Chain approximation or HNC 
closure [1], E{r) = ln(l -|- h{r) — c(r)) — h{r) -\- c[r) is the Percus-Yevick (PY) 
approximation [1]. For hard spheres, Labi'k and Malijevsky [5] introduced a 
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semiempirical approximation (LM) of E. It reproduces Monte Carlo experi- 
ments excellently. There are other approximations as the very successful clo- 
sure of Martynov and Sarkisov (MS) [6], where E[r) = \Jl -\- 2[h[r) — c(r)) — 
h[r) -\- c(r) — 1. For detailed formulas of the closures as we used them in our 
programs see below. 

Together with the closure, the OZ equation is a non-linear integral equation, 
which can be solved in general only numerically. For hard spheres in PY ap- 
proximation there is an analytic solution, too (see [7,8]). 



2 The Direct Iteration Algorithm 



The easiest algorithm for solving the OZ equation with a given closure is direct 
iteration using fast Fourier transformation. Due to the convolution theorem 
we have the following equation in k-space: 

h{k) = c{k) + p-c{k)-h{k), k = \k\. (13) 



The Fourier transforms are determined by the Fourier-Bessel transformation 
in the case of radially symmetric /(r): 

fik) = 4vr . / /(r)^^r^ dr, fir) = ^ • / fikf-^^k' dk . (14) 
J kr ZTT^ J kr 





Introducing F[r) := /(r) • r and F[k) = k ■ f{k) for f = c^h one gets the 
Fourier sine transformation 

F{k) = 47r • y" F{r) sm{kr) dr, F{r) = — ■ J F{k) sin{kr) dk . (15) 





Multiplying equation (13) by k^ and introducing V := H — C one obtains 

t = J^. (16) 

k- pC ' 



The closures can be written also in terms of c(r) = C{r)/r^ considered as 
a functional of 7(r) := r(r)/r, and the Mayer function / (see previous 
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section): 



HNC: C(r) = r ■ c(r) = r ■ {f{r) + 1) • e^MA _r(r) - r , 



PY: 


C{r) 


= ifir) + 1) • ( 


r + r(r)) — r(r) — r 


= /(r)-(r + r(r)), 


LM: 


C{r) 


= r.(/(r) + i; 


) . grCOA+i^LMC'-) -r(r) 


- , 


MS: 


C{r) 


= r.(/(r) + i; 




) — r . 



(17) 



-Elm!?") is the bridge function of Labi'k and Malijevsky [5]. 

In the case of classical one-component plasmas (see previous section) we have 
to use a somewhat different equation from Eq. (16) because of the long-range 
potential involved. Following the method of Ng [3] we obtain 



T(k) = ^ —^-C^'Hk). (18) 

k- p[c(-)ik)-km{k)) 



Here, u^^^ is the Fourier transform of the long-range part of the pair potential as 
defined in the previous section. C'^^"^ is the Fourier transform of the short-range 
part of the direct correlation function multiplied by k. The explicit relation 
to r is dependent on the closure. Here, we use only the HNC closure without 
any bridge function which is known to yield fairly good results in the region 
of plasma parameters Tp used here [3]. Then, C^*-'(r) is given by 

HNC: C(^)(r) = r exp (-/3u(^)(r) + T{r)/r) - T{r) - r (19) 



where vM"^ is the short-range part of the pair potential as defined in the previous 
section. Here, F is given by F = H — C^^\ so that the pair distribution function 
is g{r) = (F(r) + C(^)(r))/r + 1. 

Equation (16) together with any particular closure defined in Fq. (17) define 
certain integral equations. Also, equations (18) and (19) together define a fur- 
ther integral equation. The solution of any of these equations can be considered 
as a fixed point problem for the unknown function F. The integral equations 
are solved on a grid of equidistant points. Then, we put := F{i ■ AA;), 
:= F{i ■ Ar), Ar • AA; = vr/Af, where M is the number of points desired 
for calculating the former integrals. Equation (15) for the Fourier sine trans- 
formation and its inversion becomes [9] 

M-l ^ M-1 ^ 

F, = 47rAr • ^ • sin(zj-) , = — • ^ • sin(zj-) . (20) 
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Therefore, we can establish the following algorithm: Choose a r(°)(=0 for 
example) for a set of equidistant r and insert it in the closure getting C^°^. 
This can be transformed by (20) and inserted in (16) or (18) getting a T^^\ 
from which one gets r(^) by the inversion formula. This F^^^ can be used as 
a new input in the iteration process. This is done until self consistency is 
achieved, i.e. until for a given convergence threshold rj > we have 

M 

C:=J2{r^^-Tt'^) (21) 

8 = 1 

The time consuming steps are the transformations of C and F, so that it is 
desirable to reduce the number of required iterations. There are usually 200 
to 1000 iterations performed until ( < 10~^°. Therefore the aim is to use an 
acceleration method for the vector sequence F^-'^ = {Ti \ \ • • • , ^m ) 
discretized function F. 



3 Vector Extrapolation for Fixed Point Iterations 

The iterative solution of systems of nonlinear equations like the OZ equation 
(10) with some closure of the form of Eq. (12) can often be regarded as fixed 
point problems 

X = *(X) (22) 

with a parameter vector X. In the OZ case, this vector corresponds to F as 
discussed in the previous section. Such fixed point problems are often solved 
via direct iteration (Picard iteration) 

Xo , Xi = *(Xo) , . . . , = *(X„) , . . . . (23) 

In this way, a sequence of vectors X„ is generated. This sequence may or may 
not converge, and if it converges, it may or may not converge sufficiently fast. 

Especially for slowly convergent iteration sequences X„, one would like to be 
able to accelerate the convergence by some mathematical algorithms. Fortu- 
nately, this is possible. One may use vector extrapolation algorithms. These 
algorithms are a rapidly expanding field of mathematics. A good introduction 
to it is given in Chapter 4 of the textbook by Brezinski and Redivo Zaglia 
[10]. 

We discuss some general features of the acceleration of slowly convergent se- 
quences {sn}- Here, the sequence elements s„ can be numbers, vectors, ma- 
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trices et cetera. The basic principle is to use structural information hidden in 
the data. Once one has identified this structural information, it can be used 
to compute the limit faster. Usually, the result is a sequence transformation 

{sn}^=o- original sequence, (24) 
{tn}^-Q- transformed sequence. 

The transformed sequence converges hopefully in a faster way. 

The problems are to find a way to identify the type of structural information, 
and further, to construct the sequence transformation from this information. 
In order to discuss these problems, we introduce the notion of a remainder rn 
defined by 

Sn = s + rn , s = lim s„ . (25) 

Both problems are usually treated together by using a model sequence ap- 
proach. There, one takes models for the remainder r„. Then, one seeks trans- 
formations which allow - for the resulting model sequences - the exact calcu- 
lation of the limit. 

Thus, in this approach, one considers model sequences {cr„} of the form 

T 

an = cr + mn{ct,Pi) cr = r„(cr„, . . . , cr„+fc jp^) (26) 

exact 

Here, the model m„ depends on a finite number of coefficients q, and on 
further parameters pi. The transformation T eliminates the coefficients q and 
allows to calculate exactly the limit a of the model sequence {cr„} as function 
of some finite number of sequence elements cr„_|_j. The transformation T is 
specific for the model and depends parametric on the pi. 

The transformation T can also be applied to the problem sequence s„. Then, 
a sequence transformation is obtained: 

tn = Tn{sn, . . . , Sn+k\Pz) (approximate) . (27) 

The expectation is that the transformed sequence converges faster than 
the original sequence for problems that are in some sense close to the 
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model: 

s„, cr„, . (281 



A useful more special model is to factor the model remainder m„ into a re- 
mainder estimate c<;„ 7^ and a correction factor /i„(ci,7ri) according to 

Cr„ = Cr + iOnfiniCt, TTt) . (29) 

The parameters then are the tTj-. The remainder estimates cOn can also be re- 
garded as parameters. However, it is often useful to allow that the cOn depend 
also on the problem sequence. This is for instance the case for Levin's remain- 
der estimates [11] for which we display the following variants: 

''t variant": cj-n = Z\5-n_i = — -s-n—i : 

(30) 

"u variant": cOn = {n -\- l)As„_i = [n -\- l)(s„ — s„_i) . 

Here, and in the following, A denotes a difference operator acting on n in the 
form 

A/„ = - U ■ (31) 



Both t and u variants can be shown to be - up to some constant factor - good 
estimates of the true remainder r„ for large classes of sequences. In this way, 
models allow to make use of structural information. 

A further successful approach for the construction of sequence transformations 
is to use some simple basic sequence transformation To iteratively: 

To To To To 

Sn > y y . . . y ^ (32) 



This concept proved to be very successful in the case of scalar sequences [12- 
16]. 

In the vector case, one may take as the basic transformation To the transfor- 
mation 

^'^-^'^+^-"'^+MfAO„,,AOJ • ^^^^ 
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Here and in the following, quantities in capital letters like 0„, 0^^^ 
denote vectors, and (., .) denotes the usual scalar product of two vectors. The 
basic transformation (33) is exact for model sequences of the form 

E„ = E + cO„ (34) 

depending on arbitrary constants c. To apply (33) iteratively, one may take 
the same remainder estimates 0„ in each iteration. However, it is much better 
[15,16] to calculate also new remainder estimates 0^ after each iteration step 
in a hierarchical consistent [15] way. Then, one obtains the vector J transfor- 
mation as an important example that will be used in the sequel. It is a special 
case of the matrix J7 transformation that was introduced in [17]. The vector 
J7 transformation is defined by the recursive scheme 

0(0) _ a . Q(o) - Q • 




The (5^^) are numbers. They correspond to the parameters Wi in Eq. (29). Since 
they are numbers, the terms in curly braces in the recursive scheme (35) are 
also scalar numbers. 

The vector J7 transformation is a straightforward generalization of the scalar 
J7 transformation introduced in [14] that was studied intensively in [15,16,18]. 

We note that the vector J7 transformation is closely related to the vector E 
algorithm and some projection methods (See [10]). The J7 transformations 
are based on iteration. On the other hand, also remainder estimates are used 
for the J7 transformation, and in the scalar case, it is known for which model 
sequences it is exact. Hence, the J7 transformations combine features of both 
the general schemes described above. 

The problem now is how to combine vector extrapolation with Picard iteration 
for the solution of the fixed point equation (22). We describe a general method 
called cycling (cf. [10, p. 316]). It may be applied not only for the vector J7 
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transformation but for some general vector extrapolation algorithm T. In this 
method, short Picard sequences are calculated from some starting value. Then 
extrapolation is used to calculate a new starting value for the next Picard 
sequence, and so on, until convergence is achieved. This is described in more 
detail in the following. 

Thus, one constructs a vector sequence Sn from some starting vector Y by Pi- 
card iterations (cf. Eq. (23)) and uses extrapolation to compute a new starting 
vector Y' in the following way: 

50 = Y 

51 = ^{So) 

; ; (36) 

S^ = ^{Sk-i) 

Y' = T{So, Si, . . . , Sm) 



This will be called an m-cycle. In this way, direct iteration is interspersed 
with extrapolation steps that provide (hopefully better) starting values for 
the direct iteration. Also, a number d of direct iterations is normally done 
before cycling starts. Hence, the generated sequence of approximations is 

Yo, Fi 
So = Yd, Si 
So = Yd+i, Si 

So = Yd+^_i, Si = *(S'o), . . . Yd+^ = T{So, . . . , S^) cycle v 



^(Yo),... Yd = ^iYd-i) direct 
'^{So), . . . Yd+i = T{So, . . . , S^) cycle 1 
^f{So),...Yd+2 = T{So,...,S^) cycle 2 



(37) 



Thus, in each cycle a new direct iteration is started form the extrapolation 
result of the previous cycle. After the performance of cycle no. the iteration 
function ^ has been called N = d-\-v m times, and the extrapolation algorithm 
T has been applied v times. Convergence checks as described in the previous 
section can be performed at the beginning of each cycle. As noted above, the 
cycling method is applied in our case by using the vector J transformation as 
transformation T . 

In practice, one sometimes observes that for a bad starting value, there is no 
convergence of iteration methods. This holds also for Newton- Raphson-type 
approaches. Actually, in studies of deterministic models for chaotic systems 
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based on nonlinear iteration functions such a behavior is found as a generic 
case. Also, one should keep in mind that nonlinear iteration functions can 
exhibit rather peculiar behavior like fractal boundaries of the basins of attrac- 
tion or strange attractors. A further remark is that the choice of the nonlinear 
iteration function determines whether the fixed point is stable or unstable. 
For instance, there are cases where Picard iterations do not converge, but 
Newton- Raphson methods do since then a different iteration function is used. 
How does this relate to vector extrapolation? 

A first observation is that in principle, nonlinear extrapolation algorithms 
can also give rise to chaotic phenomena when the output is used iteratively 
as new input. In many applications, however, this is not the usual mode of 
operation and hence, it seems that such behavior has not been reported in the 
literature. We remark that for fixed point iterations, the combination of the 
direct iteration function ^ with cycling of a vector extrapolation method can 
be considered as a new iteration function. For instance, in the case of a 2-cycle, 
this function is vl/'(y) = T (F, '^{'^{Y))). Hence, the existence of chaotic 

phenomena for this new function is to be expected if it is nonlinear. This 
is normally the case if the old iteration function or the vector extrapolation 
algorithm are nonlinear. However, it seems that for the new iteration function, 
chaotic behavior is less probable. Put another way, it is expected that an 
unstable fixed point of ^ can become a stable fixed point of An example 
is given below. 

This is related to a second, and more important feature of extrapolation al- 
gorithms. They can transform divergent sequences into convergent ones. For 
instance, it is well-known that scalar algorithms can be used to construct ana- 
lytic continuations for power series outside of their circle of convergence in the 
form of rational approximations that converge rapidly in many cases. Also, 
suitable algorithms are able to sum divergent series as in the case of the accu- 
rate calculation of ground state energies of anharmonic oscillators from their 
strongly divergent perturbation series [19]. 

The question arises whether this can be observed also for vector extrapolation 
algorithms. As will be shown in the next section, this really is the case. Thus, 
vector extrapolation processes offer the chance to achieve convergence even 
for cases where direct iteration does not converge. 



4 Numerical Results 

Several examples have been studied, as discussed below. We used a program 
directit for direct iterations without vector extrapolation, and a program m2vj 
implementing direct iteration in combination with the u variant (see (30)) of 
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the vector J transformation defined in (35) with 



f 



(for all A;) . 



(38) 



(n + l)(n + 2) 



In the latter case, we used various numbers of direct iterations J = noffset 
as offset before cycling was started, and the length m of the m-cycles was 
varied (cf. previous section). We note that use of m2vj with m = 0, i.e., 
without cycling, means that direct iteration without extrapolation is used. 
The difference between the resulting vectors F with and without extrapolation 
is given by the quantity 



where M denotes the number of points in the grid. 

As a first example we consider a system of hard spheres at various densities 
using the PY and HNC approximation. The calculations were performed on 
a Sun Sparc workstation using 512 or 128 points for the function F at Ar = 
O.Olcr where a is the diameter of the hard sphere, respectively. The number 
of iterations that are needed to reach the convergence threshold rj (cf. Eq. 
(21)) starting from F^^^ = 0, is denoted by N. The results in Table 1 show, 
that it is possible to reduce the total number of iterations N approximately 
by a factor of two. The acceleration effect does not depend significantly on 
the choice of the approximation for the bridge function E. There was also a 
variety of values of m and noffset that yielded rather good results. As shown 
in Table 1, there is no significant difference between the resulting vectors F 
obtained by direct iteration with or without acceleration, since 6 10~^°. 

The same closures were also applied to Fennard- Jones systems. The results 
are presented in Table 2. Again, substantial reductions of the total number of 
iterations are possible. 

We also studied the Fabi'k and Malijevsky (FM) approximation and the Marty- 
nov-Sarkisov (MS) closure for hard sphere systems and for Fennard- Jones 
particles. The results are presented in Table 3. For this example, we performed 
the calculations on an Iris Indigo of Silicon Graphics using 1024 points for the 
function F, Ar = 0.005(7, and the starting vector was F^^^ = 0. As noted above, 
the case m = corresponds to performing only direct iterations without any 
acceleration. 

As a further example, we present results of calculations on classical one- 
component plasmas using Ng renormalization [3] with HNC closure (without 




(39) 
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Table 1 

Hard spheres 



Diameter: a = 1, number density: p* : = 


pla-\ 


Starting 


vector: r(°) 


= 0, total 


num- 


ber of iterations: iV, ^ 


^rid size: 


M, threshold: rj 


= 1-10- 


-1° (cf. Eq. ( 


21)), deviation: 


S (cf. Eq. (39)), cycle 


length: 


m, offset: 


'ioffset, calculation on Sun workstation. 


case program 


m 


^^offset 










N 


A directit 


- 












592 


m2vj 


8 


14 


no convergence up to 900 iterations 


m2vj 


9 


14 


7.092 


1 n — 10 
• 10 






295 


m2vj 


9 


15 


4.955 


1 n-io 
• 10 






306 


m2vj 


10 


14 


4.yo5 


• lU 






4UU 


m2vj 


10 


15 


o.z4U 


1 n-io 
• lU 






/I /? 

o4d 


B directit 


- 












1 O /I 

124 


m2vj 


7 


10 


4.yU4 


• lU 






51 


m2vj 


8 


5 


5.640 


• 10-11 






41 


m2vj 


8 


10 


4:.U01 


• 10-1° 






47 


C directit 


- 












488 


m2vj 


8 


10 


2.318 


•10-9 






209 


D directit 




— 










384 


m2vj 


9 


14 


3.080 


•10-9 






455 


m2vj 


10 


13 


2.716 


•10-9 






179 


m2vj 


10 


14 


2.963 


•10-9 






257 


m2vj 


11 


13 


3.022 


•10-9 






182 


A: PY^ approximation, p* = 


= 0.7 (Ar 


= 0.01 


(7, M = 


512) 






B: PY^ approximation, p* = 


: 0.4 (Ar 


= 0.01 


(7, M = 


512) 






C: PY^ approximation, p* = 


: 0.7 (Ar 


= 0.04 


(7, M = 


128) 






D: HNC^ approximation, p* 


= 0.7 (Ar = 0.01(7, M 


= 512) 






^ The abbreviations 


for the closures . 


are explained in Section 1 


and Eq. i 


(17). 



taking any further bridge diagrams into account). The results are displayed in 
Table 4. Again, we used M = 1024 points in the grid and started from T^^^ = 0. 
It is worth noting that for higher plasma parameters Tp it is advantageous to 
choose a higher number noffset of iterations without any acceleration than at 
low Tp, but the cycle length is not to be changed very much. Actually we 
found our best results nearly at all values of Tp at m = 10. This is favorable 
because the computational costs for the extrapolation rise with increasing m. 
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Table 2 

Lennard-Jones potential 



LJ+ parameters: ci = 1, /3e 
Ar = 0.01(7, threshold: r] = 
workstation, other symbols 


= 0.5, number density: p* = p/a~'^. Grid: M 
1 • 10-1°, starting vector: r(°) = 0, calculation 
see Table 1. 


= 512, 
on Sun 


case 


program 


m 


^^offset 


S 




N 


A 


directit 










137 




fn2vj 


7 


10 


6.532 • 


10-1° 


51 




m2vj 


8 


10 


4.950 • 


10-1° 


56 


B 


directit 










9^8 




m2vj 


7 


20 


8.174 • 


10-9 


349 




m2vj 


7 


30 


4.296 • 


10-9 


383 




m2vj 


8 


20 


4.705 • 


10-9 


372 




m2vj 


10 


50 


7.854 • 


10-9 


480 


C 


directit 










615 




m2vj 


8 


30 


1.234 • 


10-1° 


517 



A: PY^ approximation, p* = 0.5 
B: PY^ approximation, p* = 0.9 
C: HNC^ approximation, p* = 0.9 

^ The abbreviations for the closures are explained in Section 1 and Eq. (17). 
+ The Lennard-Jones potential is defined in Eq. (5). 



As one can see from the tables, the acceleration is in any case successful 
if one performs a sufficient number of iterations (usually approximately 10) 
without acceleration before cycling. At higher number densities as in Table 
3 it can be necessary to perform more direct iterations. Here, the best result 
was obtained after 100 or more direct iterations. But see in contrast the last 
example of Table 3 at a very high number density, where the best result is 
obtained using no iterations without acceleration. The length of cycles m is not 
needed to be very high, usually 8 to 20. This is desirable because in this way 
storage and computing times for the cycling are reduced. The use of longer 
cycles can even be disadvantageous as one can see from Tables 2 or 3. Eor 
Lennard-Jones systems the results are of similar quality as in the hard sphere 
case (see Tables 2 and 3). Especially, note the relatively high number densities 
up to p* = 1.2 in Table 3, which was considered up to now to be intractable 
with direct iteration [20], when starting from r(°) = 0. 
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Table 3 



Hard 


spheres and Lennard-Jones 


potential 






Sphere diameter/LJ+ parameter: a = 


1. Starting 


vector: 


: 1^") = 0, threshold: r] = 


1 • 10" 


■1° (cf. Eq. (21)), : 


grid size: M 


= 1024. Program: 


■m2vj, calculation on Iris 


Indigo 


, other symbols see Table 1. 








case 




m 




^^offset 


N 


A 











576 






14 




150 


361 






15 




100 


309 






16 




120 


308 






16 




150 


389 


B 











913 






9 







611 






9 




1 


572 






15 




100 


581 






15 




200 


537 


C 











840 






15 




20 


661 






15 




100 


469 






15 




120 


681 


D 











956 






7 







241 






7 




1 


378 



A: LM^ bridge function, hard spheres, p* = 0.75 (Ar = 0.005(t) 

B: MS^ approximation, hard spheres, p* = 0.80 (Ar = 0.005(t) 

C: MS^ approximation, Lennard-Jones+, /3e = 0.5, p* = 0.90 (Ar = 0.005(t) 

D: MS^ approximation, Lennard-Jones+, /3e = 0.1, p* = 1.20 (Ar = O.Olcr) 

^ The abbreviations for the closures are explained in Section 1 and Eq. (17). 

+ The Lennard-Jones potential is defined in Eq. (5). 

In order to assess the additional costs for the extrapolation steps in the cy- 
cling algorithm, we measured the total CPU time to run our programs for 
several examples. As it is well-known, such measurements have to be inter- 
preted cautiously since the results depend not only on the basic algorithm, 
but also on the skills of the programmer and the machine architecture and 
utilization when the programs are run. 
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Table 4 

One-component Plasmas 

Plasma parameter: Fp, a = 1 (cf. Eq. (7)). HNC^^ approximation. Ng renormaliza- 
tion, a (parameter in error function): 1.08/a (cf. Eq. (8) and [3]). Number density: 
p = 3/(47ra3). Grid: M = 1024, Ar = O.Ola. Starting vector: r(°) = 0. Program: 
■m2vj, calculation on Iris Indigo. Other symbols see Table 1. 





m 


^^offset 


N 


10 







70 


10 


5 





55 


10 


7 


5 


30 


10 


8 


5 


33 


10 


10 





23 


50 







252 


50 


7 


10 


75 


50 


10 





78 


50 


10 


10 


77 


100 







458 


100 


8 


40 


275 


100 


8 


60 


412 


100 


10 


10 


275 


1 nn 

lUU 


1 n 

lU 


ou 


101 


100 


14 


40 


341 


120 







537 


120 


8 


100 


614 


120 


10 


50 


414 


120 


10 


100 


189 


120 


12 


100 


361 


^ The abbreviations 


for the closures 


are explained in Section 1 and Eq. 


(17). 



The basic result is that for typical cycle lengths and grid sizes the costs per 
iteration for the extrapolation part are of the same order as the costs for the 
direct iteration alone. Incidentally, this shows that for the direct iteration the 
costs are rather low. This is due to the use of the fast Eourier transform. Eor 
more complicated iteration functions, the relative costs of the extrapolation 
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are expected to be better. But even for the very fast direct iterations of the 
current apphcation, there are cases where it is possible to reduce the total 
execution time by a factor of up to two by using vector extrapolation. This 
follows from the following examples. 

The case of classical one-component plasmas is treated in Table 5. For given 
values of the plasma parameter Fp, the total time taken by direct iteration 
only, i.e., with m = 0, is compared to the total time for the case with the lowest 
number of total iterations in Table 4) for m 7^ 0, i.e., including acceleration. 
For example, for Tp = 120, the CPU time for m = is compared to the CPU 
time for m = 10 and noffset = 100. 

Table 5 
CPU Times 



Plasma parameter (cf. Eq. (7)): Tp. Calculations on Sun workstation 





Direct iteration alone 


Direct iteration + acceleration 


Saving 


100 


37.7 s 


28.0 s 


25.7 % 


120 


44.1 s 


22.0 s 


50.1 % 



A similar result was obtained for the last example of Table 3: 88.5 seconds 
of the direct iteration versus 44.6 seconds with acceleration, a saving of 49.6 
percent of the CPU time. 



We want to stress that - for the proposed method - the computing times are 
very low. On workstations, each run can be completed in times of the order of 
a minute. This should be compared to computer simulations where computing 
times of the order of hours or more are required. 

Finally, we want to give an example that shows that there are cases where 
vector extrapolation can be used successfully to find fixed points even when 
the direct iteration does not converge when started from F^^^ = 0. Actually, 
for this starting value, even the Newton- Raphson-type algorithm of Fabi'k, 
Malijevsky, and Vonka [20] does not converge. For the latter, we use a Fortran 
program called lensub. 

The example is a hard-sphere system for p* = 0.85 with LM closure. Using a 
grid of size M = 512 and with Ar = O.Olcr, the results are the following: 

- The Newton- Raphson-type algorithm [20] does not converge for a starting 
vector F(°) = 0. It does converge for better starting vectors (for instance 
from a run with lower density). 

- The direct iteration does not converge for a starting vector F^^^ = 0. It 
also does not converge when the solution obtained with lensub is used as 
starting vector. This reveals that this starting vector corresponds to an 
unstable fixed point for the direct iteration. See also Figure 1 and 2. 
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- The direct iteration in combination with the vector extrapolation converges 
lor a starting vector F^^^ = to the same solution as the Newton- Raphson- 
type algorithm. Some examples are displayed in Table 6. 

Table 6 

Hard Spheres 

Density p* = 0.85. LM"!" closure. Program: m2fj, calculation on Iris Indigo. Grid: 
M = 512, Ar = 0.01(7. Threshold: rj = 10"^°. Other symbols see Table 1. 



Example m f^offset ^ 

1 9 40 991 

2 9 100 891 

3 9 150 881 

4 9 200 981 



^ The abbreviations for the closures are explained in Section 1 and Eq. (17). 

In Eigure 1, we plot semi-logarithmically the values of ( defined in Eq. (21) for 
this example for the direct iteration. Thus, ( measures the distances between 
consecutive vectors in the iteration. 



1 1 1 1 r 




Fig. 1. Unstable Fixed Point of the Direct Iteration 

It is clearly seen that the direct iteration first seems to converge but then it 
is starting to change rapidly until a quasiperiodic behavior is reached. The 
latter is displayed in an expanded representation in Eigure 2. 

In Eigure 3, we plot the values of ( defined in Eq. (21) for the case of Example 
3 in Table 6, i.e., for the direct iteration in combination with vector extrapo- 
lation. Convergence is rather smooth apart from some step-like features. 
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InC 



T 1 1 1 1 1 1 1 r 



3.5 



2.5 



1.5 



900 905 910 915 920 925 930 935 940 945 950 

N 

Fig. 2. Quasiperiodic Behavior of the Direct Iteration 




200 



300 



400 



500 

N 



600 



700 



800 



900 



Fig. 3. Convergence of Vector-extrapolated Iteration 

Before leaving this topic we note that for the same system and grid for a 
somewhat smaller density of p* = 0.8 the direct iteration converges from a 
starting vector F^^^ = while the Newton- Raphson-type algorithm of [20] 
does not converge using this starting vector. 



An important advantage of the present algorithm is the very simple form 
of the direct iteration. This makes the treatment of more complicated cases 
as particles with dipoles or more realistic model potentials defined on the 
grid very easy. Essentially, one only has to reprogram the subroutine for the 
computation of the Mayer function. In the case of a renormalization of the 
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OZ scheme as for OCPs, the changes are also easy to implement. All other 
algorithms as that due to Labi'k, Malijevsky, and Vonka [20, p. 710] need 
more sophisticated programming, than the method presented here, because 
no calculation of the gradient of the iteration function is required as in all 
algorithms based on the Newton- Raphson method. 

The promising results of this work suggest to study the combination of direct 
iteration in combination with acceleration methods also for more complicated 
model potentials and multicomponent systems. The applicability of further 
known vector extrapolation processes in the field of the Ornstein-Zernike equa- 
tion should be studied. In the opinion of the authors, also the development of 
new powerful vector extrapolation processes is possible and desired that are 
tailored to speed up vector iteration processes even further. Thus, we stress 
that the methods of the present work still have the potential to further im- 
provements. Finally, it is an interesting question whether it is possible to use 
acceleration methods profitably also in combination with more complicated 
algorithms like Newton- Raphson iterations. In summary, in the context of 
the Ornstein-Zernike equation, the combination of direct iteration methods 
with vector extrapolation has been shown to be a fruitful alternative to other 
methods. 
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